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Percolation Approach to Study Connectivity in 
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Abstract. We study neural connectivity in cultures of rat hippocampal neurons. We measure the 
neurons' response to an electric stimulation for gradual lower connectivity, and characterize the size 
of the giant cluster in the network. The connectivity undergoes a percolation transition described 
by the critical exponent j3 ~ 0.65. We use a theoretic approach based on bond-percolation on a 
graph to describe the process of disintegration of the network and extract its statistical properties. 
Together with numerical simulations we show that the connectivity in the neural culture is local, 
characterized by a gaussian degree distribution and not a power law one. 
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^ INTRODUCTION 

^: 

Neurons in living networks form a highly rich web of connections in which activity 
_ flows between neurons through synapses. The most fascinating living neural network is 

\Q • the human brain, but its complex architecture, functionality, and computation capability 

is still far from being fully understood. More impressive is that the 100 billion neurons 
are not randomly connected, but rather form elaborate circuits with specific tasks. Con- 
nectivity thus appears as the fundamental feature to understand the potential of a living 
neural network. 



o 

O , Unravelling the detailed connectivity diagram of a living neural network is a painstak- 

ing process. For a brain, a small section of it, or even for a small neural culture, with 
~ 10^ neurons and ~ 10^ connections in just 1 mm^, this task is, at present, unfeasible. 

KA \ In the brain, substantial progress has been attained in the description of the connectiv- 

J-j \ ity in the mammalian cortex [1, 2], or the analysis of brain functional networks [3, 4]. 

However, only in the small invertebrate C. elegans [5] it has been possible to map out, 
in a Herculean project, the connectivity of its 302 neurons. It is not surprising then that 
other approaches, different than the pure physiological ones, are being introduced to ex- 
tract information about the connectivity of neural networks or, at least, some relevant 
statistical properties. 

Biological neural networks have caught the attention of Physicists and Mathemati- 
cians following the "burst" of interest that complex networks and random graphs have 
experienced in the last decade [6, 7]. Graph theory has permitted to reduce the com- 
plexity of a rich variety of natural and artificial networks (e.g. Internet, e-mail, social, 
collaborations, or genetic networks) in terms of basic concepts that retain their most 
important features, such as the presence of a power law connectivity, clustering, or the 
small world phenomena. One of these concepts, which is related with percolation theory 
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FIGURE 1. (a) Phase contrast image of a small region of a neural culture. Spherical objects are neurons. 
(b) Fluorescence image. Bright spots are cell bodies. The scale bar is 50 jUm in both images, (c) Sketch of 
the experimental setup. The F{t) plot shows an example of the fluoresce signal of a spiking neuron as a 
function of time. The vertical dashed Une indicates the excitation time. 



[8, 9], is the characterization of the giant cluster (or giant component) of the network 
and how it disintegrates as links or nodes are removed. A power law connectivity for 
instance makes the network robust to random attacks, but vulnerable to directed attacks, 
since the removal of just a small number of highly connected nodes destroys the gi- 
ant component [9]. The problem of resilience is of great interest for biological neural 
networks, and makes the study of neural connectivity of enormous importance. 

Next, we will see how concepts of graph and percolation theory can be used to extract 
statistical information about the connectivity in living neural networks. We will describe 
our experimental results on connectivity in neural cultures and their analysis in terms of 
bond-percolation on a graph [10]. Together with numerical simulations of the model we 
show that the connectivity in neural cultures is characterized by a Gaussian distribution 
(and not a power law one), and with the presence of some locality and clustering. 

EXPERIMENTAL SETUP AND PROCEDURE 



Experiments (see Ref. [10] for details) were performed on primary cultures of rat 
hippocampal neurons, that are plated on glass coverslips (Fig. la). Neurons develop 
dendrites and axons shortly after plating, creating a dense web of connections in a few 
days (Fig. lb). Cultures were used 14-20 days after plating, when the network is fully 
developed and its activity is governed by the balance between excitatory and inhibitory 
neurons. About 20% of the neurons are known to be inhibitory [11]. 

Neurons were electrically stimulated through bath electrodes (Fig. Ic), and the cor- 
responding voltage drop V measured with an oscilloscope. Neuronal activity was moni- 
tored using fluorescence calcium imaging, and data processed to record the fluorescence 
intensity F as a function of time. Neural spiking activity is detected as a sharp increase 
of the fluorescence intensity. 

The connectivity of the network was gradually weakened by blocking the AMPA glu- 
tamate receptors of excitatory neurons with the receptor antagonist CNQX. We studied 
the role of inhibition by either leaving active or blocking the GABA receptors with the 
corresponding antagonist bicuculine. For simplicity, we label the network containing 
both excitatory and inhibitory neurons by Gei, and the network with excitatory neurons 
only by Ge- The response of the network for a given CNQX concentration was quan- 
tified as the fraction of neurons $ that fired in response to the electric stimulation at 



voltage V. Response curves 4>(y) were obtained by increasing the stimulation voltage 
from 2 to 6 V in steps of 0.1 — 0.5 V. At the end of the experiments, the culture was 
washed of CNQX to verify that the initial network connectivity was recovered. 

MODEL 

We consider a simplified model of the network in terms of bond-percolation on a graph. 
The neural network is represented by the directed graph G. Our main simplifying as- 
sumption is the following: A neuron has a probability / = f{V) to fire as a direct 
response to the electric excitation, and it always fires if any one of its input neurons fire 
(Fig. 2a). This approach ignores the fact that more than one input is needed to excite 
a neuron, and that connections are gradually weakened rather than abruptly removed. 
However, the aim of the model is to provide the simplest scenario to understand the ex- 
perimental observations, and not the actual, highly complex behavior of neural cultures. 
/ is the natural unit in which to measure the response of the network, and by a change 
of variable the measured response curves 4>(y) can be expressed as $(/). 

The fraction of neurons in the network that fire for a given value of / defines the firing 
probability 4>(/). 4>(/) increases with the connectivity of G, because any neuron along 
a directed path of inputs may fire and excite all the neurons downstream (Fig. 2a). All 
the upstream neurons that can thus excite a certain neuron define its input-cluster or 
excitation-basin. It is therefore convenient to express the firing probability as the sum 
over the probabilities ps of a neuron to have an input-cluster of size s— I (Figs. 2b-c), 

$(/) = / + ( 1 — f)P (any input neuron fires) 

oo oo 

= /+(i-/)iAv(i-(i-/rM = i-iAs.(i-/r, (1) 

5=1 ^ ^ .S=l 

where we used the probability conservation Y^sPs = 1- ^(/) increases monotonically 
with / and ranges between $(0) = and 4>(1) = 1. The deviation of 4>(/) from linearity 
manifests the connectivity of the network (for disconnected neurons 4>(/) = /). Eq. (1) 
indicates that the observed firing probability 4>(/) is actually one minus the generating 
function H{x) (or the z-transform) of the cluster-size probability ps [12], 

oo 

H{x) = Y,Psx' = l-^{f), (2) 

i=i 

where x = 1 — /. One can extract from H{x) the input-cluster size probabilities ps, 
formally by the inverse z-transform, or more practically, in the experiments, by fitting 
H{x) to a polynomial in x. 

Once a giant component emerges (Fig. 2d) the observed firing pattern is significantly 
altered. In an infinite network, the giant component always fires no matter what the firing 
probability / is. This is because even a very small / is sufficient to excite one of the 
infinitely many neurons that belong to the giant component. We account for this effect 
by splitting the neuron population into a fraction g that belongs to the giant component 
and always fires and the remaining fraction 1 — ^ that belongs to finite clusters. This 
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FIGURE 2. (a) Percolation model. The neuron represented in grey fires either in response to an external 
excitation or if any of its input neurons fire. At the highest connectivity, this neuron has input-clusters 
s — 1 = (self-excitation), 7 (left branch), 6 (right branch), and 13 (both branches). At lower connectivity, 
its input-clusters are reduced to sizes and 3. (b) Corresponding Ps{s) distributions, obtained by counting 
all input-clusters for all neurons. Insets: H{x) functions for the Ps{s) distributions (solid lines), compared 
with independent neurons, H{x) ~ x (dashed lines), (c) Example of the sensitivity of ps{s) to loops. Left: 
neurons forming a chain-like connectivity give a Ps{s) distributed uniformly. Center: closing the loop 
by adding just one link collapses Ps{s) to a single peak. Right: additional links increase the average 
connectivity {k), but do not modify /7i(.s). (d) Concept of giant component. The grey areas outline the size 
of the giant component g (biggest cluster) for gradually smaller connectivity c. 



modifies the summation on cluster sizes into 

oo 

^(/) =^ + (1-^) [/+(1 -/)P(any inp. neu. fires)] = 1 -(1 -^) £p,(l -/)\ (3) 

5=1 

As expected, at the limit of almost no self-excitation / — t- only the giant component 
fires, $(0) = g, and $(/) monotonically increases to 4>(1) = 1. With a giant component 
present the relation between H{x) and the firing probability changes, obtaining 



H{x) = Y,Psx' 



1-j? 



(4) 



The size of the giant component decreases with the connectivity. At a critical connectiv- 
ity Co the giant component disintegrates and its size is comparable to the average cluster 
size in the neural network. This behavior corresponds to a percolation transition, sep- 
arating a system of small, fragmented clusters to one with a fast growing giant cluster 
that comprises most of the network. 

EXPERIMENTAL RESULTS 



Examples of the response curves ^{V) for Gei and Ge networks are shown in Figs. 3a 
and 3b. At one extreme, with [CNQX] = the network is fully connected. All neurons 
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FIGURE 3. (a) and (b) Examples of response curves 4>(y) for Gei (top) and Ge (bottom) networks at 
different concentrations of CNQX. Grey bars indicate the size of the giant component. Lines are a guide 
to the eye except for 1 /xM and 10 jiM that are fits to error functions, (c) Size of the giant component as 
a function of the connectivity c for Gei networks (circles) and Ge networks (squares). Lines are a guide 
to the eye. Some CNQX concentrations are indicated for clarity, (d) Log-log plot of the power law fits 
g - |1 -c/q/, with Co = 0.36 ±0.02, j8 = 0.66 ±0.05 for Gei , and Co = 0.24 ±0.02, j3 = 0.63 ±0.05 
for Ge ■ (e) H{x) functions for the response curves shown in (a) and for [CNQX] > 100 nM. Lines are 
polynomial fits up to order 20. (f) Corresponding cluster size distribution Ps{s)- 



form a single cluster that comprises the entire network. A few neurons with low firing 
threshold suffice to activate the entire culture, leading to a very sharp response curve. 
At the other extreme, with high concentrations of CNQX (~ 10 /iM) the network is 
completely disconnected, and the response curve is given by the individual neurons' 
response. 4>(y) for individual neurons (denoted as ^^{V)) is well described by an error 

function 4>(V) = 0.5 + 0.5 erf f 7- " j . This indicates that the firing threshold of a neuron 

in the network follows a gaussian distribution with mean Vq and width 2oq. 

Intermediate CNQX concentrations induce partial blocking of the synapses. Some 
neurons break off into separated clusters, while a giant cluster still contains most of 
the remaining neurons. The response curves are then characterized by a big jump that 
corresponds to the biggest cluster {giant component), and two tails that correspond 
to smaller clusters of neurons with low or high firing threshold. Beyond a critical 
concentration (around 500 nM for Gei networks and 700 nM for Ge networks) a giant 
component cannot be identified and the whole response curve is then also well described 
by an error function. 

The biggest cluster in the network characterizes the giant component g. For each 
response curve, g is measured as the biggest fraction of neurons that fire together in 
response to the electric excitation, as shown by the grey bars in Figs. 3a and 3b. The 
size of the giant component was studied as a function of the connectivity probability (or 



synaptic strength) between two neurons [10], given by c = 1/(1 + [CNQX]/^^), with 
Kfi = 300 nM, and takes values between (full blocking) and 1 (full connectivity). 

The breakdown of the network for both Gei and Ge networks is shown in Fig. 3c. The 
giant component for Gei networks breaks down at much lower CNQX concentrations 
compared with Ge networks, indicating that the effect of inhibition on the network is 
to effectively reduce the number of inputs that a neuron receives on average. The be- 
havior of the giant component indicates that the neural network undergoes a percolation 
transition, described by the power law g ~ 1 1 — c/co\^ . Power law fits for Gei and Ge 
networks give the same j8 ^ 0.65 within the experimental error (Fig. 3d), indicating that 
j8 is an intrinsic property of the network. 

Finally, we have studied the size distribution Ps{s) for clusters that do not belong to 
the giant component. Ps{s) has been obtained by constructing the experimental function 
H{x) and after fitting a polynomial Y.sPs^- Since / = ^oo(y) is the response curve 
for individual neurons (Figs. 3a and 3b) and x = 1 — /, the function H{x) for each 
response curve is obtained by plotting 1 — 4>(y) as a function of 1 — 4>oo(V). For curves 
with a giant component present, its contribution is eliminated and the resulting curve 
normalized by the factor 1 — ^. Fig. 3e shows the H{x) functions for the response curves 
of Fig. 3a. The corresponding Ps{s) distribution, obtained from fits up to order 20, is 
shown in Fig. 3f. Overall, the clusters start out relatively big to rapidly become smaller 
for gradually higher concentrations of CNQX. Ps{s) is characterized by isolated peaks, 
indicating that loops and strong locality may be present in the neural culture. An example 
that illustrates the strong effect of loops on Ps{s) is shown in Fig. 2c. Since Ps{s) is 
obtained by fitting polynomials on H{x), the accuracy in the description of Ps{s) is 
limited by the resolution ofH{x) which, in turn, is limited by the experimental resolution 
in 4>(y). In addition, since Ps{s) is a probability distribution, the fit is carried out with 
two constraints, reducing the freedom of fitting: the p^ coefficients have to be positive 
and their sum has to be one. Hence, the Ps{s) distribution presented in Fig. 3f shows the 
correct behavior, but not the precise details of the distribution of input-clusters. 

NUMERICAL SIMULATIONS 

The model has been derived from classic bond percolation theory and has an analytic 
solution that yields precise results. However, the model contains a series of simplifying 
assumptions that may have an effect on the results. The numerical simulations that 
we present next are oriented to investigate the effect of removing or relaxing these 
assumptions, and to provide a physical picture for the connectivity in the network. 

Three assumptions of the model are unrealistic. First, it assumes that one input 
suffices to activate a neuron, while in reality a number of input neurons must spike 
for the target neuron to fire. Second, the effect of CNQX is to bind and block AMPA 
glutamate receptor molecules, and consequently to continuously reduce the synaptic 
strength, so that bonds are in reality gradually weakened rather than abruptly removed. 
Third, the model assumes a tree-like connectivity, while in the living culture loops and 
clusters may exits. The numerical simulations have been applied to test that none of 
these assumptions change the main results of the model, i.e. that the giant component 
undergoes a percolation transition at a critical connectivity cq, and that the analysis of 



H{x) provides the distribution of input-clusters in the network. 

The numerical simulations also provide the framework to study different degree 
distributions and their effect in the critical exponent /3. A Gaussian distribution gives 
j8 ^ 0.66, as in the experiments, while a power law distribution, pk{k) ~ k^ , gives j8 
equal to or larger than one, where its exact value depends on the exponent A [13]. 



Numerical method 

The neural network was simulated as a directed random graph G{N,kj/Q) in which 
each vertex is a neuron and each edge is a synaptic connection between two neurons 
[14]. The graph was generated by assigning to each edge an input/output connectivity 
kj/o according to a predetermined degree distribution. Next, a connectivity matrix Qy 
was generated by randomly connecting pairs of neurons with a link of initial weight 
1 until each vertex was connected to kj/Q links. The process of gradual weakening of 
the network was simulated in one case by removing edges, and in the second case by 
gradually reducing the bond strength from 1 to 0. The connectivity c is defined for the 
case of removing bonds as the fraction of remaining edges, and for the case of weakening 
bonds as the bond strength. 

Each neuron has a threshold v, to fire in response to the external voltage, and all neu- 
rons have a threshold T to fire in response to the integrated input from their neighbors. 
Since the experiments show that the probability distribution for independent neurons to 
fire in response to an external voltage is Gaussian, the V;'s are distributed accordingly. 
For the simple case of removing links, the global threshold T differentiates networks 
where a single input suffices to excite a target neuron from those where multiple inputs 
are necessary. When links are weakened T plays a more subtle role, and determines the 
variable number of input neurons that are necessary to make a target neuron spike. 

The state of each neuron, inactive (0) or active (1) was kept in a state vector S. In the 
first simulation step, a neuron fires in response to the external voltage if the "excitation 
voltage" V is greater than its individual threshold v,, i.e. V > v, — )■ 5,- = 1. 

In the subsequent simulation steps, a neuron fires due to the internal voltage if the 
integration over all its inputs at a given iteration is larger than T: Y^CijSj >T^Si=\. 
The simulation iterates until no new neurons fire. The network response 4>(y) is then 
measured as the fraction of neurons that fired during the stimulation. The process is 
repeated for increasing values of V , until the entire network gets activated, ^{V) = 1. 
Then, the network is weakened and the exploration in voltages started again. 



Simulations results 

Analysis of the model 

To study the validity of the model, we have first considered 4 different situations: 
removing or weakening edges, and for T = 1 or T = 5. In all cases the connectivity 
is set to be Gaussian for both input and output degree distributions. The results of the 
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FIGURE 4. Numerical simulations for 4 different cases. Shown are the response curves ^(V), the cor- 
responding H{x) functions (inset), and the characterization of the percolation transition for (a) removing 
edges, T=l; (b) weakening edges, T=l; (c) removing edges, T=5; and (d) weakening edges, T=5. 



simulation are presented in Fig. 4. All 4 studied cases give qualitatively similar results, 
with response curves ^{V) that are comparable to the ones observed experimentally, and 
with a giant component clearly identifiable. The analysis of the percolation transition 
gives /3 ~ 0.66 in all 4 cases, in agreement with the value measured experimentally. As 
expected, the simulations with weakening bonds and for T = 5 (five spiking neurons 
required to excite the target neuron) provide the response curves that are more similar to 
the ones observed experimentally. However, it is remarkable that the simplest case of the 
model (breaking bonds with T = 1) already gives valid results. This indicates that, with 
the limitations of the model, the percolation approach proves to be remarkably powerful 
in describing the behavior observed experimentally. 

The other important assumption of the model is the effect of the presence of loops 
in the network. Although loops are very rare in a random (Erdos-Renyi) graph, the 
connectivity in neural cultures is not random, and locality and neighboring probably 
may play an important role. However, graph theory tells us that most loops will be found 
in the giant component, where all neurons anyway light up and their effect is therefore 
irrelevant to our analysis. Clusters outside the giant component are in general tree-like, 
and thus the important analysis to be considered is what happens when finite clusters do 
have loops. The simulations show that the response curves and the percolation transition 
are not significantly altered if loops are allowed, providing similar results to the ones 
shown in Fig. 4. 

Loops, however, do affect the clusters size distribution Ps{s), which is then character- 
ized by the presence of isolated peaks. To explore to which extent Ps{s) was sensitive 
to loops, we performed simulations considering different levels of clustering. The first 
graph that we analyzed was one with an artificially induced highly clustered connec- 



tivity. We generated a network where most of the links are located in highly connected 
clusters, with only weak connections between clusters. The Ps{s) distribution obtained 
from the breakdown of the connectivity in such a clustered network showed that the 
position of the dominant peaks corresponded to the size of the highly connected clus- 
ters. Next, having demonstrated the importance of highly connected clusters, we went 
on to consider realizations of the graph that would be more similar to the experimental 
network. To do that, we introduced the notion of geometry and of distance, placing all 
vertices on a spatial grid. Three different configurations were used: (i) a Gaussian con- 
nectivity with no locality, (ii) a Gaussian connectivity with local connections and (iii) a 
Gaussian connectivity with distance dependent link strength. For the first case no domi- 
nant peaks where identifiable. For the second one the existence of isolated peaks is more 
apparent. But for the third case the reinforced connectivity significantly increases the 
probability to have isolated input-clusters, similar to what we observe experimentally. 



Role of inhibition and analysis ofH(x) 

We have studied the role of inhibition in the network by randomly selecting a sub- 
group of nodes and assigning them negative weights to simulate inhibitory neurons. 
Then, simulations with the same conditions described above were repeated and different 
excitation/inhibition ratios explored. The results indicate that the critical exponent /3 is 
independent of the balance between excitation and inhibition, in agreement with the ex- 
perimental observations. The results also show that the critical connectivity cq at which 
the giant component disintegrates does depend on the number of inhibitory neurons, and 
that this is a linear dependence. 

Finally, we have verified with the simulations that the cluster distribution Ps{s) ob- 
tained from the polynomial fit of H{x) does not differ significantly from the Ps{s) dis- 
tribution directly extracted from the connectivity matrix Ctj. Small deviations are a con- 
sequence of the constraints Y.Ps = ^ and < Ps < 1 in the polynomial fits, and in the 
uncertainty in removing the contribution of the giant component in the H{x) functions. 
This analysis gives validity to the Ps{s) distribution measured experimentally. 

DISCUSSION AND CONCLUSIONS 

By comparing the exponent j8 measured experimentally with the one obtained from 
the simulations we conclude that the connectivity in the neural culture is Gaussian. 
Simulations, however, are based on a random graph, while the real neural network 
is not, and one may think that the neural culture is actually better described by a 
two-dimensional, lattice-like network. Percolation on two-dimensional lattices gives 
a critical exponent /3 ~ 0.14, independent on the lattice structure. The value of the 
exponent increases rapidly with the dimensionality of the lattice, with j8 ~ 0.41 and 0.64 
for three and four dimensions, respectively. In a system described by a 2-D structure, 
additional dimensions can be viewed as a gradual increase of long-range correlations. 

The physical picture that we think may exist in the neural culture is that neurons 
are essentially connected to their neighbors, but with some long-range correlations. 



Axons can easily extend 300 jim in a neural culture, connecting neurons as far as 30 
cell bodies. The concept that locality is important is in fact quite natural when one 
thinks of the nature of the culture. Neurons are distributed homogeneously over the 
glass, and most likely all neurons start to form connections at the same time and at 
the same rate. This hints at a structure where neurons are highly connected with their 
neighbors. This is also suggested by the distribution of input-clusters Ps{s), which shows 
that neurons are highly connected between them even after the giant component has 
begun disintegrating, forming local clusters with a significant presence of loops. We 
have also seen that neurons surrounded by many others tend to fire first in response to 
the external excitation, and that aggregates of neurons tend to fire together, with their 
collective response maintained even when the connectivity is reduced. 

In summary, we have presented experimental results on the connectivity in neural 
cultures, and showed that connectivity undergoes a percolation transition characterized 
by a critical exponent /3 ^ 0.65. The experimental results were studied in the framework 
of percolation on a graph, and extracted the distribution of connected components in the 
network. Numerical simulations of the model were used to construct a physical picture of 
the connectivity in the neural network, and showed that the connectivity is characterized 
by a Gaussian degree distribution, with strong locality and clusterization. 
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